bayesian hw 2

chapter 5

Author

Congyuan He

Published

November 28, 2024

EXERCISE 4.15.

Find Jeffreys’ prior for \(\theta\) based on a random sample of size \(n\) when (a) \(y_i|\theta ∼ Pois(\theta)\), (b) \(y_i|\theta ∼ Exp(\theta)\), (c) \(y_i|\theta ∼ Weib(2,\theta)\), (d) \(y_i|\theta\)is negative binomial as in Example 4.3.3.

Answer:

Definte

In general, for a one parameter problem, Fisher’s information is defined to be the expected value of the negative of the second derivative of the log-likelihood. Jeffreys’ prior is defined as being proportional to the square root of the Fisher information.

参数化不变性：Jeffreys’ prior 的设计保证了在重新参数化（如从换成 = g()）时，先验分布形式保持一致。这种不变性使其成为真正的非信息性先验。
非信息性：它试图提供一种客观的先验选择，当缺乏明确的先验知识时，使用 Jeffreys’ prior 不会引入人为的偏倚。
依赖模型：Jeffreys’ prior 依赖于数据生成模型的结构，因为 Fisher 信息量取决于似然函数。

(a) \(y_i|\theta ∼ Pois(\theta)\)

We have, \[ f(y|\theta) = \theta^y e^{-\theta }/y! \]

then,

\[\begin{align} L(\theta|y) &= \frac{\theta^{\sum y} e^{-\theta y}}{\prod y!} \\ log[L(y |\theta)] &= \sum y * ln\theta - \theta *\sum y - ln(\prod y_i !) \\ log[L(y |\theta)]' &= \frac{\sum y}{\theta} + \sum y \\ log[L(y |\theta)]'' &= - \frac{\sum y}{\theta^2} \\ -E[log[L(y |\theta)]''] &= -E[- \frac{\sum y}{\theta^2} ] \\ J(\theta) &= \frac{n \theta }{\theta^2} \\ f(\theta) &\propto \theta^{-0.5} \end{align}\]

(b) \(y_i|\theta ∼ Exp(\theta)\)

we have,

\[ f(y|\theta) = \theta e^{−θy}I(0,\infty)(y) \]

then,

\[\begin{align} L(y |\theta) &= \prod \theta e^{-\theta y} \\ log[L(y |\theta)] &= nln\theta - \theta \sum y \\ log[L(y |\theta)]' &= \frac{n}{\theta} - \sum y \\ log[L(y |\theta)]'' &= - \frac{n}{\theta^2} \\ -E[log[L(y |\theta)]''] &= -E[- \frac{n}{\theta^2}] \\ J(\theta) &= \frac{n}{\theta^2} \\ f(\theta) &\propto \theta^{-1} \end{align}\]

(c) \(y_i|\theta ∼ Weib(2,\theta)\)

we have,

\[ f(y|\alpha,\theta) = 2 \theta y * exp (−\theta y^{2}) I(0,∞)(y) \]

then,

\[\begin{align} L(y |\theta) &= 2\theta ^n \prod y exp(-n \theta \sum y^2) \\ log[L(y |\theta)] &= nln2\theta + ln\prod y -n\theta \sum y^2 \\ log[L(y |\theta)]' &= \frac{2n}{2\theta} -n \sum y^2 \\ log[L(y |\theta)]'' &= - \frac{n}{\theta^2} \\ -E[log[Ln(y |\theta)]''] &= -E[- \frac{n}{\theta^2}] \\ J(\theta) &= \frac{n}{\theta^2} \\ f(\theta) &\propto \theta^{-1} \end{align}\]

(d) \(y_i|\theta\) is negative binomial

we have,

\[ L(\theta｜y_2) = \binom{y_2-1}{s-1}\theta^{s}(1-\theta)^{y_2-s} \]

then,

\[\begin{align} Log(L(\theta｜y_2)) &\propto sln(\theta) + (y_2 -s)ln(1-\theta) \\ Log(L(\theta｜y_2))' &\propto \frac{s}{\theta}-\frac{y_2-s}{1-\theta} \\ Log(L(\theta｜y_2))'' &\propto - \frac{s}{\theta^2}-\frac{y_2-s}{(1-\theta)^2} \\ -E[Log(L(\theta｜y_2))''] &\propto -E[- \frac{s}{\theta^2}-\frac{y_2-s}{(1-\theta)^2}] \\ J(\theta) &\propto \frac{y_2}{\theta} + \frac{y_2}{1-\theta}\\ p(\theta)&\propto \sqrt{\frac{1}{\theta(1-\theta)}} \end{align}\]

手写修改后的答案：

EXERCISE 5.2.

suppressMessages(library(rjags))
library(coda)

model_string <- "
  model {
    y ~ dbin(theta, n)
    theta ~ dbeta(12.05, 116.06)
  }
"

data_list <- list(n = 2430, y = 219)

jags_model <- jags.model(textConnection(model_string), 
                         data = data_list, 
                         n.chains = 1,    
                         n.adapt = 5000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 1
   Unobserved stochastic nodes: 1
   Total graph size: 5

Initializing model

update(jags_model, n.iter = 5000)

samples <- coda.samples(jags_model, 
                        variable.names = c("theta"), 
                        n.iter = 5000)  

summary(samples)


Iterations = 10001:15000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean             SD       Naive SE Time-series SE 
     9.043e-02      5.781e-03      8.176e-05      1.047e-04 

2. Quantiles for each variable:

   2.5%     25%     50%     75%   97.5% 
0.07956 0.08642 0.09028 0.09423 0.10202

# plot(samples)

# Plot the Beta(12.05, 116.06) distribution
plot(seq(0, 1, length.out = 1000), dbeta(seq(0, 1, length.out = 1000), 12.05, 116.06), 
     type = "l", col = "#4169E1", lwd = 2, 
     main = "Beta(12.05, 116.06) vs Beta(231.05, 2327.06)", 
     ylim = c(0,70),
     xlab = expression(theta), ylab = "Density")

# Add the Beta(231.05, 2327.06) distribution
lines(seq(0, 1, length.out = 1000), dbeta(seq(0, 1, length.out = 1000),231.05, 2327.06), 
      col = "red", lwd = 2, lty = 2)  # Use dashed line type for differentiation

# Add a legend
legend("topright", legend = c("Beta(12.05, 116.06)", "Beta(231.05, 2327.06)"), 
       col = c("#4169E1", "red"), lwd = 2, lty = c(1, 2))

EXERCISE 5.3.

Using calculus, find the mode and 5th percentile of a Beta(10,1) distribution.

\[\begin{align} beta(10,1) &= 10\theta^9 \\ f(\theta)'&=90\theta^8 \end{align}\] 单调递增 mode = 1

\[\begin{align} q&=\int_{0}^{c} f(\theta)d\theta \\ q&= c^{10} \end{align}\]

5 分位数为 \(0.05^{1/10}\)

EXERCISE 5.4.

Using calculus, find a and b such that a Beta(a, b) distribution has a mode of 1 and a 5th percentile of 0.2. mode=\((a-1)/(a+b-2)\),令其=1，则b=1

EXERCISE 5.5.

Derive formula (1), including the formula for θ0.

\[\begin{align} \theta_0 &= \frac{a-1}{a+b-2} \\ a(\theta_0 -1) &= 2\theta_0 - b\theta_0 -1 \\ a &= \frac{2\theta_0 - b\theta_0 - 1}{\theta_0 -1} \end{align}\]

EXERCISE 5.6.

Use BetaBuster to find the Beta(a, b) priors for mode 0.75 and 5th percentile 0.60, and for mode 0.01 and 99th percentile 0.02. What is the Beta prior when the mode is 1 and the first percentile is 0.80?

EXERCISE 5.7.

The distributions θ ∼ Beta(1.6, 1) and θ ∼ Beta(1, 0.577) both have a mode of 1. Find Pr[θ < 0.5] analytically for each. Does BetaBuster give the appropriate parameters for the Beta distributions?

pbeta(0.5,1.6,1)

[1] 0.329877

pbeta(0.5,1,0.577)

[1] 0.3296437

EXERCISE 5.8.

model_string <- "
  model{ 
  gamma[1] ~ dbeta(a1,b1) 
  gamma[2] ~ dbeta(a2,b2) 
  theta <- w*gamma[1] + (1-w)*gamma[2]
  w ~ dbern(p) }
"

data_list <- list(a1 = 10, b1 = 20, a2 = 20, b2 = 10, p=0.5 )

jags_model <- jags.model(textConnection(model_string), 
                         data = data_list, 
                         n.chains = 1,    
                         n.adapt = 5000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 0
   Unobserved stochastic nodes: 3
   Total graph size: 13

Initializing model

update(jags_model, n.iter = 5000)

samples <- coda.samples(jags_model, 
                        variable.names = c("theta"), 
                        n.iter = 5000)  

# summary(samples)

hist(samples[[1]])

model_string <- "
  model{ 
  gamma[1] ~ dbeta(a1,b1) 
  gamma[2] ~ dbeta(a2,b2) 
  theta <- w*gamma[1] + (1-w)*gamma[2]
  w ~ dbern(p) }
"

data_list <- list(a1 = 10, b1 = 20, a2 = 20, b2 = 10, p=0.2 )

jags_model <- jags.model(textConnection(model_string), 
                         data = data_list, 
                         n.chains = 1,    
                         n.adapt = 5000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 0
   Unobserved stochastic nodes: 3
   Total graph size: 13

Initializing model

update(jags_model, n.iter = 5000)

samples <- coda.samples(jags_model, 
                        variable.names = c("theta"), 
                        n.iter = 5000)  

# summary(samples)

hist(samples[[1]])

model_string <- "
  model{ 
  gamma[1] ~ dbeta(a1,b1) 
  gamma[2] ~ dbeta(a2,b2) 
  theta <- w*gamma[1] + (1-w)*gamma[2]
  w ~ dbern(p) }
"

data_list <- list(a1 = 1, b1 = 1, a2 = 20, b2 = 10, p=0.5 )

jags_model <- jags.model(textConnection(model_string), 
                         data = data_list, 
                         n.chains = 1,    
                         n.adapt = 5000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 0
   Unobserved stochastic nodes: 3
   Total graph size: 13

Initializing model

update(jags_model, n.iter = 5000)

samples <- coda.samples(jags_model, 
                        variable.names = c("theta"), 
                        n.iter = 5000)  

# summary(samples)

hist(samples[[1]])

model_string <- "
  model{ 
  gamma[1] ~ dbeta(a1,b1) 
  gamma[2] ~ dbeta(a2,b2) 
  theta <- w*gamma[1] + (1-w)*gamma[2]
  w ~ dbern(p) }
"

data_list <- list(a1 = 0.1, b1 =0.2, a2 = 20, b2 = 10, p=0.5 )

jags_model <- jags.model(textConnection(model_string), 
                         data = data_list, 
                         n.chains = 1,    
                         n.adapt = 5000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 0
   Unobserved stochastic nodes: 3
   Total graph size: 13

Initializing model

update(jags_model, n.iter = 5000)

samples <- coda.samples(jags_model, 
                        variable.names = c("theta"), 
                        n.iter = 5000)  

# summary(samples)

hist(samples[[1]])

EXERCISE 5.9.

各种调整参数测试

model_string <- "
model{
theta ~ dbeta(a,b)T(,t)
}
"
data_list <- list(a=2,b=3,t=0.9)
inits <- list(theta = 0.5)
jags_model <- jags.model(textConnection(model_string),
                      data=data_list,
                      n.chains=1,
                      n.adapt = 1000,
                      inits = inits)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 0
   Unobserved stochastic nodes: 1
   Total graph size: 4

Initializing model

update(jags_model,n.iter = 1000)
samples <- coda.samples(jags_model,
                        variable.names = c('theta'),
                        n.iter = 5000)
hist(samples[[1]])

EXERCISE 5.12.

model_string <- "
model{
y[1] ~ dbin(theta[1],n[1]) 
y[2] ~ dbin(theta[2],n[2]) 

# betabuster
theta[1] ~ dbeta(3.2846,6.3307) 

theta[2] ~ dbeta(1.5317,1.5317) 


odds[1] <- theta[1]/(1-theta[1]) 
odds[2] <- theta[2]/(1-theta[2]) 
RD <- theta[2]-theta[1] 
RR <- theta[2]/theta[1] 
OR <- odds[2]/odds[1] 
test <- step(RD)
}
"

data_list <- list(n=c(134,400), y=c(54,224)) 

jags_model <- jags.model(textConnection(model_string),
                      data=data_list,
                      n.chains=1,
                      n.adapt = 5000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 2
   Unobserved stochastic nodes: 2
   Total graph size: 18

Initializing model

update(jags_model,n.iter = 5000)

samples <- coda.samples(jags_model,
                        variable.names = c('theta','RD','RR','OR','test'),
                        n.iter = 5000)
summary(samples)


Iterations = 10001:15000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

           Mean      SD  Naive SE Time-series SE
OR       1.9586 0.39258 0.0055519      0.0068229
RD       0.1606 0.04761 0.0006733      0.0008332
RR       1.4175 0.16106 0.0022777      0.0028135
test     0.9996 0.02000 0.0002828      0.0002828
theta[1] 0.3990 0.04080 0.0005769      0.0007176
theta[2] 0.5595 0.02449 0.0003463      0.0004350

2. Quantiles for each variable:

            2.5%    25%    50%    75%  97.5%
OR       1.29120 1.6842 1.9240 2.1915 2.8174
RD       0.06363 0.1291 0.1616 0.1926 0.2514
RR       1.13587 1.3052 1.4081 1.5138 1.7629
test     1.00000 1.0000 1.0000 1.0000 1.0000
theta[1] 0.32121 0.3715 0.3975 0.4261 0.4798
theta[2] 0.51106 0.5431 0.5594 0.5758 0.6066

对比发现区别不大，受数据影响大

EXERCISE 5.14.

model_string <- "
model{
y[1] ~ dbin(theta[1],n[1]) 
y[2] ~ dbin(theta[2],n[2]) 

theta[1] ~ dunif(0,1) 
theta[2] ~ dunif(0,1) 

odds[1] <- theta[1]/(1-theta[1]) 
odds[2] <- theta[2]/(1-theta[2]) 

OR <- odds[1]/odds[2] 
}
"

data_list <- list(n=c(7,16), y=c(6,8)) 

jags_model <- jags.model(textConnection(model_string),
                      data=data_list,
                      n.chains=1,
                      n.adapt = 5000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 2
   Unobserved stochastic nodes: 2
   Total graph size: 13

Initializing model

update(jags_model,n.iter = 5000)

samples <- coda.samples(jags_model,
                        variable.names = c('OR'),
                        n.iter = 5000)
summary(samples)


Iterations = 10001:15000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean             SD       Naive SE Time-series SE 
        7.5689        15.6947         0.2220         0.2452 

2. Quantiles for each variable:

   2.5%     25%     50%     75%   97.5% 
 0.6377  2.0738  3.9610  7.7736 36.4600

# 用共轭结论
# 均匀分布beta(1,1)
nsim <- 10000
a <- b <- 1
n1 <- 7; y1 <- 6; n2 <- 16; y2 <- 8

theta1_post <- rbeta(nsim,a+y1,b+n1-y1)
theta2_post <- rbeta(nsim,a+y2,b+n2-y2)

OR <- theta1_post/(1-theta1_post)/(theta2_post/(1-theta2_post))
summary(OR)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
   0.1189    2.1084    4.0621    8.1318    8.0350 1035.2793

关联并不能代表因果

EXERCISE 5.15.

model_string <- "
model{ 
for(i in 1:2){ y[i] ~ dbin(theta[i],n[i]) } 
theta[2] ~ dbeta(a,b) 
delta ~ dnorm(mu, prec) 
theta[1] <- exp(delta)*theta[2]/(1-theta[2]*(1-exp(delta))) 
OR <- theta[1]/(1-theta[1])/(theta[2]/(1-theta[2])) }
"

data_list <- list(n=c(7,16), y=c(6,8),a=0.5,b=0.5,mu=2, prec=1/2)

jags_model <- jags.model(textConnection(model_string),
                      data=data_list,
                      n.chains=1,
                      n.adapt = 10000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 2
   Unobserved stochastic nodes: 2
   Total graph size: 22

Initializing model

update(jags_model,n.iter = 10000)

samples <- coda.samples(jags_model,
                        variable.names = c('OR'),
                        n.iter = 10000)
summary(samples)


Iterations = 20001:30000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 10000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean             SD       Naive SE Time-series SE 
       12.7363        19.3095         0.1931         0.3205 

2. Quantiles for each variable:

  2.5%    25%    50%    75%  97.5% 
 1.270  3.913  7.351 14.303 57.329

mean(samples[[1]]>1)

[1] 0.9872

mean(samples[[1]]>2)

[1] 0.927

# 换先验 dunif(log(0.02), log(50))
model_string <- "
model{ 
for(i in 1:2){ y[i] ~ dbin(theta[i],n[i]) } 
theta[2] ~ dbeta(a,b) 
delta ~ dunif(log(0.02), log(50)) # 这里改了
theta[1] <- exp(delta)*theta[2]/(1-theta[2]*(1-exp(delta))) 
OR <- theta[1]/(1-theta[1])/(theta[2]/(1-theta[2])) }
"

data_list <- list(n=c(7,16), y=c(6,8),a=1,b=1)

jags_model <- jags.model(textConnection(model_string),
                      data=data_list,
                      n.chains=1,
                      n.adapt = 10000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 2
   Unobserved stochastic nodes: 2
   Total graph size: 24

Initializing model

update(jags_model,n.iter = 10000)

samples <- coda.samples(jags_model,
                        variable.names = c('OR'),
                        n.iter = 10000)
summary(samples)


Iterations = 20001:30000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 10000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean             SD       Naive SE Time-series SE 
       11.0003        10.9282         0.1093         0.1493 

2. Quantiles for each variable:

   2.5%     25%     50%     75%   97.5% 
 0.9082  3.2367  6.9098 14.8292 42.4630

mean(samples[[1]]>1)

[1] 0.9681

mean(samples[[1]]>2)

[1] 0.8728

EXERCISE 5.16.

model_string <- "
model {
  for (i in 1:2) {
    y[i] ~ dbin(theta[i], n[i]) 
  }

  theta_t1 ~ dbeta(13.3221, 6.2809)  
  theta_t2 ~ dbeta(6.2809,13.3221)              
  gamma ~ dbeta(2, 2)                

  theta[1] <- theta_t1 * gamma / (theta_t1 * gamma + theta_t2 * (1 - gamma))  
  theta[2] <- (1 - theta_t1) * gamma / ((1 - theta_t1) * gamma + (1 - theta_t2) * (1 - gamma))  
  
  OR <- theta[1]/(1-theta[1])/(theta[2]/(1-theta[2])) }

"

data_list <- list(n=c(7,16), y=c(6,8))

jags_model <- jags.model(textConnection(model_string),
                      data=data_list,
                      n.chains=1,
                      n.adapt = 10000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 2
   Unobserved stochastic nodes: 3
   Total graph size: 27

Initializing model

update(jags_model,n.iter = 10000)

samples <- coda.samples(jags_model,
                        variable.names = c('OR'),
                        n.iter = 10000)
summary(samples)


Iterations = 20001:30000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 10000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean             SD       Naive SE Time-series SE 
       5.92988        4.04559        0.04046        0.05042 

2. Quantiles for each variable:

  2.5%    25%    50%    75%  97.5% 
 1.598  3.319  4.877  7.301 16.424

mean(samples[[1]]>1)

[1] 0.9972

mean(samples[[1]]>2)

[1] 0.9442

model_string <- "
model {
  for (i in 1:2) {
    y[i] ~ dbin(theta[i], n[i]) 
  }

  theta_t1 ~ dbeta(28.3393, 42.009)  
  theta_t2 ~ dbeta(42.009,28.3393)              
  gamma ~ dbeta(2, 2)                

  theta[1] <- theta_t1 * gamma / (theta_t1 * gamma + theta_t2 * (1 - gamma))  
  theta[2] <- (1 - theta_t1) * gamma / ((1 - theta_t1) * gamma + (1 - theta_t2) * (1 - gamma))  
  
  OR <- theta[1]/(1-theta[1])/(theta[2]/(1-theta[2])) }

"

data_list <- list(n=c(7,16), y=c(6,8))

jags_model <- jags.model(textConnection(model_string),
                      data=data_list,
                      n.chains=1,
                      n.adapt = 10000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 2
   Unobserved stochastic nodes: 3
   Total graph size: 27

Initializing model

update(jags_model,n.iter = 10000)

samples <- coda.samples(jags_model,
                        variable.names = c('OR'),
                        n.iter = 10000)
summary(samples)


Iterations = 20001:30000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 10000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean             SD       Naive SE Time-series SE 
      0.620732       0.201921       0.002019       0.002476 

2. Quantiles for each variable:

  2.5%    25%    50%    75%  97.5% 
0.3131 0.4779 0.5910 0.7305 1.0977

mean(samples[[1]]>1)

[1] 0.0494

mean(samples[[1]]>2)

[1] 0

EXERCISE 5.17.

没有患病数据做不出来

# 验证5.16的先验对不对
a1 <- 13.32;b1 <- 6.28
a2 <- 6.28;b2 <- 13.32
# mode
(a1-1)/(a1+b1-2);(a2-1)/(a2+b2-2)

[1] 0.7

[1] 0.3

qbeta(0.05,a1,b1)

[1] 0.499983

qbeta(0.95,a2,b2)

[1] 0.500017

nism <- 10000
theta_tilde1 <- rbeta(nism,a1,b1)
theta_tilde2 <- rbeta(nism,a2,b2)
gamma <- rbeta(nsim,2,2)

theta1 <- theta_tilde1 * gamma/(theta_tilde1 * gamma+theta_tilde2 *(1-gamma))
theta2 <- (1 - theta_tilde1) * gamma / ((1 - theta_tilde1) * gamma + (1 - theta_tilde2) * (1 - gamma)) 

hist(theta1)

hist(theta2)

plot(theta1,theta2)

EXERCISE 5.18.

直接把两个数据相加，用原始的先验

# 重复5.14 计算出 theta1 2 的后验分布
model_string <- "
model{
y[1] ~ dbin(theta[1],n[1]) 
y[2] ~ dbin(theta[2],n[2]) 

theta[1] ~ dunif(0,1) 
theta[2] ~ dunif(0,1) 

}
"

data_list <- list(n=c(7,16), y=c(6,8)) 

jags_model <- jags.model(textConnection(model_string),
                      data=data_list,
                      n.chains=1,
                      n.adapt = 5000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 2
   Unobserved stochastic nodes: 2
   Total graph size: 8

Initializing model

update(jags_model,n.iter = 5000)

samples <- coda.samples(jags_model,
                        variable.names = c('theta[1]','theta[2]'),
                        n.iter = 5000)
summary(samples)


Iterations = 10001:15000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

           Mean     SD Naive SE Time-series SE
theta[1] 0.7743 0.1343 0.001900       0.002837
theta[2] 0.4982 0.1173 0.001659       0.002153

2. Quantiles for each variable:

           2.5%    25%    50%    75%  97.5%
theta[1] 0.4591 0.6948 0.7982 0.8772 0.9673
theta[2] 0.2664 0.4178 0.4976 0.5813 0.7223

# 更新

# 用之前的后验当作新的先验
plot(density(samples[[1]][,1]))

plot(density(samples[[1]][,2]))

library(MASS)

# 对 samples[[1]][,1] 的密度拟合 Beta 分布
data1 <- samples[[1]][,1]
fit1 <- fitdistr(data1, densfun = "beta", start = list(shape1 = 1, shape2 = 1))

Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

# 对 samples[[1]][,2] 的密度拟合 Beta 分布
data2 <- samples[[1]][,2]
fit2 <- fitdistr(data2, densfun = "beta", start = list(shape1 = 1, shape2 = 1))

Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

# 打印拟合结果
fit1$estimate # 包含 shape1 和 shape2 的值

  shape1   shape2 
6.834040 1.994196

fit2$estimate

  shape1   shape2 
8.552017 8.616468

model_string <- "
model{
y[1] ~ dbin(theta[1],n[1]) 
y[2] ~ dbin(theta[2],n[2]) 

theta[1] ~ dbeta(7.061013, 1.999425) 
theta[2] ~ dbeta(9.053316, 8.994726) 

odds[1] <- theta[1]/(1-theta[1]) 
odds[2] <- theta[2]/(1-theta[2]) 

OR <- odds[2]/odds[1] 
}
"

data_list <- list(n=c(37,75), y=c(35,32)) 

jags_model <- jags.model(textConnection(model_string),
                      data=data_list,
                      n.chains=1,
                      n.adapt = 5000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 2
   Unobserved stochastic nodes: 2
   Total graph size: 16

Initializing model

update(jags_model,n.iter = 5000)

samples <- coda.samples(jags_model,
                        variable.names = c('OR'),
                        n.iter = 5000)
summary(samples)


Iterations = 10001:15000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean             SD       Naive SE Time-series SE 
     0.0793840      0.0458026      0.0006477      0.0009964 

2. Quantiles for each variable:

   2.5%     25%     50%     75%   97.5% 
0.01892 0.04646 0.07003 0.10131 0.19544

敏感性分析换不同的先验结果差不多

EXERCISE 5.23.

model_string <- "
model{ 
for(i in 1:n){ y[i] ~ dnorm(mu, tau) }  # 似然函数

# mu ~ dflat() # 会报错 不能自动产生初始值
# mu ~ dnorm(0, 1/1000000) 
mu ~ dnorm(4.75, 1/0.0163) 
tau ~ dgamma(c,d) 
sigma <- 1/sqrt(tau)
gamma <- phi((4.4-mu)/sqrt(1/tau)) 
prob <- step(4.4 - y[13]) }  
"

data_list <- list(y=c(4.20,4.36,4.11,3.96,5.63,4.50, 5.64,4.38,4.45,3.67,5.26,4.66,NA),
                  c=0.001, d=0.001,n=13) 


jags_model <- jags.model(textConnection(model_string),
                      data=data_list,
                      n.chains=1,
                      n.adapt = 5000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 12
   Unobserved stochastic nodes: 3
   Total graph size: 32

Initializing model

update(jags_model,n.iter = 5000)

samples <- coda.samples(jags_model,
                        variable.names = c('mu','sigma','y[13]','gamma','prob'),
                        n.iter = 5000)
summary(samples[[1]])


Iterations = 5001:10000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

        Mean      SD  Naive SE Time-series SE
gamma 0.3293 0.06304 0.0008915      0.0008271
mu    4.6889 0.10576 0.0014956      0.0015272
prob  0.3318 0.47091 0.0066596      0.0066596
sigma 0.6662 0.15309 0.0021651      0.0022494
y[13] 4.6844 0.70362 0.0099507      0.0099507

2. Quantiles for each variable:

        2.5%    25%    50%    75%  97.5%
gamma 0.2050 0.2864 0.3308 0.3729 0.4533
mu    4.4812 4.6187 4.6882 4.7580 4.9034
prob  0.0000 0.0000 0.0000 1.0000 1.0000
sigma 0.4414 0.5620 0.6406 0.7411 1.0333
y[13] 3.2806 4.2366 4.6786 5.1283 6.0725

plot(samples[[1]][,'gamma'])

model_string <- "
model{ 
for(i in 1:n){ y[i] ~ dnorm(mu, tau) } 
 mu ~ dnorm(0, 1/1000000) 
# mu ~ dnorm(4.75, 1/0.0163) 
tau ~ dgamma(c,d) 
sigma <- 1/sqrt(tau)
gamma <- phi((4.4-mu)/sqrt(1/tau)) 
prob <- step(4.4 - y[13]) }  
"

data_list <- list(y=c(4.20,4.36,4.11,3.96,5.63,4.50, 5.64,4.38,4.45,3.67,5.26,4.66,NA),
                  c=0.001, d=0.001,n=13) 


jags_model <- jags.model(textConnection(model_string),
                      data=data_list,
                      n.chains=1,
                      n.adapt = 5000)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 12
   Unobserved stochastic nodes: 3
   Total graph size: 32

Initializing model

update(jags_model,n.iter = 5000)

samples <- coda.samples(jags_model,
                        variable.names = c('mu','sigma','y[13]','gamma','prob'),
                        n.iter = 5000)
summary(samples[[1]])


Iterations = 5001:10000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

        Mean     SD Naive SE Time-series SE
gamma 0.4002 0.1097 0.001552       0.001470
mu    4.5702 0.1996 0.002823       0.002747
prob  0.3952 0.4889 0.006915       0.006915
sigma 0.6774 0.1567 0.002216       0.002215
y[13] 4.5890 0.7121 0.010070       0.010070

2. Quantiles for each variable:

        2.5%    25%    50%    75%  97.5%
gamma 0.1991 0.3224 0.3972 0.4741 0.6239
mu    4.1749 4.4427 4.5687 4.6999 4.9645
prob  0.0000 0.0000 0.0000 1.0000 1.0000
sigma 0.4508 0.5656 0.6536 0.7594 1.0481
y[13] 3.1776 4.1269 4.5785 5.0362 6.0551

plot(samples[[1]][,'gamma'])

EXERCISE 5.24.

Solve for b above when the lower 5th percentile is specified by the expert. \[ a − 1.645 \sqrt \frac{1}{b} = l \]

90% upper

\[ a + 1.281552 \sqrt \frac{1}{b} = u \]

90% lower \[ a - 1.281552 \sqrt \frac{1}{b} = l \]

find the bs for u = 70 and l = 58

a <- 65
b <- 1/0.0163
u <- 70
l <- 58
b_u = (1.281552/(u - a))^2

b_l = (1.281552/(a - l))^2

b_u

[1] 0.06569502

b_l

[1] 0.03351787

EXERCISE 5.25.

在设置先验的时候，\(\mu\) 单独一个分布 \(\tau\) 单独一个分布，两者是独立的
后验分布更新的时候两者概率也是独立的\(p(μ,τ) = p(μ)p(τ).\)

EXERCISE 5.26.

这里建议先看原文推导

#1 R Code for finding prior on sigma 
alpha <- 0.90 
beta <- 0.95 
a <- 65 # Best guess for mu 
tildegamma <- 85 # Best guess for gamma_alpha
tildeu <- 91 # Best guess percentile of gamma_alpha 
zalpha <- 1.28 # qnorm(0.90,0,1) 
f <- 2.588 # Initial value for f 
# Could use a sequence of values, say f <- seq(1,50,1) 
sigma0 <- (tildegamma - a)/zalpha
e <- 1 + sigma0*f 
# We must find the Gamma(e,f) distribution that
# has beta-percentile = tildesigmabeta 
tildesigmabeta <- (tildeu - a)/zalpha 
trialq <- qgamma(beta,e,f) # Return beta-percentile for the 
                           # selected gamma distribution 
trialq # If trialq = tildesigmabeta

[1] 20.31031

tildesigmabeta # stop and pick corresponding f

[1] 20.3125

[1] 41.4375

[1] 2.588

#2 R Code for finding prior on tau 
alpha <- 0.90
beta <- 0.95
zalpha <- 1.28 
a <- 65 
tildegamma <- 85 
tildel <- 79  # 最后一问这里改成70
d <- 1328
tau0 = (zalpha/(tildegamma - a))^2
c = 1 + tau0*d 
tildetaubeta = (zalpha/(tildel - a))^2 
trialq = qgamma(beta,c,d) 
trialq

[1] 0.008358841

tildetaubeta

[1] 0.008359184

[1] 6.439488

plot(seq(0, 60, length.out = 1000), 
     dgamma(seq(0, 100, length.out = 1000), shape = 41.4375, rate = 2.588), 
     type = "l", main = "Density of σ ~ Gamma(41.4375, 2.588)", col = "blue", lwd = 2)

tau_samples <- rgamma(10000, shape = 6.439, rate = 1328)
sigma_samples <- (1 / tau_samples)^0.5

# 绘制 σ 的密度图
lines(density(sigma_samples), col = "red", lwd = 2)

# 添加图例
legend("topright", legend = c(" σ ~ Gamma(41.4375, 2.588)", 
                               "τ ~ Gamma(6.439, 1328)"),
       col = c("blue", "red"), lwd = 2)

your best guess for the mean exam score is 60,

#1 R Code for finding prior on sigma 
alpha <- 0.90 
beta <- 0.95 
a <- 60 # Best guess for mu 改这里
tildegamma <- 85 # Best guess for gamma_alpha
tildeu <- 91 # Best guess percentile of gamma_alpha 
zalpha <- 1.28 # qnorm(0.90,0,1) 
f <- 3.068 # Initial value for f   # 二分法手调 懒得写函数了
# Could use a sequence of values, say f <- seq(1,50,1) 
sigma0 <- (tildegamma - a)/zalpha
e <- 1 + sigma0*f 
# We must find the Gamma(e,f) distribution that
# has beta-percentile = tildesigmabeta 
tildesigmabeta <- (tildeu - a)/zalpha 
trialq <- qgamma(beta,e,f) # Return beta-percentile for the 
                           # selected gamma distribution 
trialq # If trialq = tildesigmabeta

[1] 24.21879

tildesigmabeta # stop and pick corresponding f

[1] 24.21875

[1] 60.92188

a=60
u=70
b <- (1.645/(u-a))^2 
b

[1] 0.02706025

Note

\(\mu ～ N(60, 1/0.02706025)\)

\(\sigma ～ Gamma(60.92188, 3.068)\)

you are 95% sure that the mean exam score is less than 65,

a=60
u=65
b <- (1.645/(u-a))^2 
b

[1] 0.108241

#1 R Code for finding prior on sigma 
alpha <- 0.90 
beta <- 0.95 
a <- 60 # Best guess for mu 
tildegamma <- 85 # Best guess for gamma_alpha
tildeu <- 91 # Best guess percentile of gamma_alpha 
zalpha <- 1.28 # qnorm(0.90,0,1) 
f <- 3.068 # Initial value for f   # 二分法手调 懒得写函数了
# Could use a sequence of values, say f <- seq(1,50,1) 
sigma0 <- (tildegamma - a)/zalpha
e <- 1 + sigma0*f 
# We must find the Gamma(e,f) distribution that
# has beta-percentile = tildesigmabeta 
tildesigmabeta <- (tildeu - a)/zalpha 
trialq <- qgamma(beta,e,f) # Return beta-percentile for the 
                           # selected gamma distribution 
trialq # If trialq = tildesigmabeta

[1] 24.21879

tildesigmabeta # stop and pick corresponding f

[1] 24.21875

[1] 60.92188

Note

\(\mu ～ N(60, 1/0.108241)\)

\(\sigma ～ Gamma(60.92188, 3.068)\)

your best guess for the 90th percentile of exam scores is 80, and

a=60
u=65
b <- (1.645/(u-a))^2 
b

[1] 0.108241

#1 R Code for finding prior on sigma 
alpha <- 0.90 
beta <- 0.95 
a <- 60 # Best guess for mu 
tildegamma <- 80 # Best guess for gamma_alpha
tildeu <- 91 # Best guess percentile of gamma_alpha 
zalpha <- 1.28 # qnorm(0.90,0,1) 
f <- 0.932# Initial value for f   # 二分法手调 懒得写函数了
# Could use a sequence of values, say f <- seq(1,50,1) 
sigma0 <- (tildegamma - a)/zalpha
e <- 1 + sigma0*f 
# We must find the Gamma(e,f) distribution that
# has beta-percentile = tildesigmabeta 
tildesigmabeta <- (tildeu - a)/zalpha 
trialq <- qgamma(beta,e,f) # Return beta-percentile for the 
                           # selected gamma distribution 
trialq # If trialq = tildesigmabeta

[1] 24.21494

tildesigmabeta # stop and pick corresponding f

[1] 24.21875

[1] 15.5625

Note

\(\mu ～ N(60, 1/0.108241)\)

\(\sigma ～ Gamma(15.5625, 0.932)\)

you are 95% sure that the 90th percentile is less than 90

a=60
u=65
b <- (1.645/(u-a))^2 
b

[1] 0.108241

#1 R Code for finding prior on sigma 
alpha <- 0.90 
beta <- 0.95 
a <- 60 # Best guess for mu 
tildegamma <- 80 # Best guess for gamma_alpha
tildeu <- 90 # Best guess percentile of gamma_alpha 
zalpha <- 1.28 # qnorm(0.90,0,1) 
f <- 1.089 # Initial value for f   # 二分法手调 懒得写函数了
# Could use a sequence of values, say f <- seq(1,50,1) 
sigma0 <- (tildegamma - a)/zalpha
e <- 1 + sigma0*f 
# We must find the Gamma(e,f) distribution that
# has beta-percentile = tildesigmabeta 
tildesigmabeta <- (tildeu - a)/zalpha 
trialq <- qgamma(beta,e,f) # Return beta-percentile for the 
                           # selected gamma distribution 
trialq # If trialq = tildesigmabeta

[1] 23.43242

tildesigmabeta # stop and pick corresponding f

[1] 23.4375

[1] 18.01562

Note

\(\mu ～ N(60, 1/0.108241)\)

\(\sigma ～ Gamma(18.01562, 1.089)\)

b中对应部分修改为70

 tildel <- 79  # 最后一问这里改成70

EXERCISE 5.29.

run_jags_model <- function(n1, mu1, tau1, n2, mu2, tau2, a1, b1, c1, d1, a2, b2, c2, d2) {

  y <- rnorm(n1, mu1, sqrt(1/tau1))
  x <- rnorm(n2, mu2, sqrt(1/tau2))
  
  model_string <- "
  model {
    for (i in 1:n[1]) {
      y[i] ~ dnorm(mu[1], tau[1])
    }
    for (j in 1:n[2]) {
      x[j] ~ dnorm(mu[2], tau[2])
    }
    
    for (r in 1:2) {
      mu[r] ~ dnorm(a[r], b[r])
      tau[r] ~ dgamma(c[r], d[r])
      sigma[r] <- sqrt(1/tau[r])
    }
    
    meandiff <- mu[1] - mu[2]
    sdratio <- sigma[1] / sigma[2]
    
    prob[1] <- step(meandiff) # Pr(meandiff > 0 | data)
    prob[2] <- step(sdratio - 1) # Pr(sdratio > 1 | data)
  }"
  
    data_list <- list(y = y, x = x, n = c(n1, n2), 
                    a = c(a1, a2), b = c(b1, b2), 
                    c = c(c1, c2), d = c(d1, d2))
  
  
  # 运行JAGS模型
  jags_model <- jags.model(textConnection(model_string), data = data_list, 
                            n.chains = 1, n.adapt = 5000)
  
  update(jags_model , n.iter = 5000)
  
  # 提取结果
  results <- coda.samples(jags_model, variable.names = c("mu", "tau", "meandiff", "sdratio", "prob"), n.iter = 5000)
  
  return(results)
}

set.seed(2024) 
n1 <- 30     # 第一组的样本大小
mu1 <- 5     # 第一组的均值
tau1 <- 1    # 第一组的精度

n2 <- 20     # 第二组的样本大小
mu2 <- 3     # 第二组的均值
tau2 <- 0.5  # 第二组的精度

# 两个分布的参考先验参数
a1 <- 0; b1 <- 0.001; c1 <- 0.001; d1 <- 0.001
a2 <- 0; b2 <- 0.001; c2 <- 0.001; d2 <- 0.001

# 运行JAGS模型
results <- run_jags_model(n1, mu1, tau1, n2, mu2, tau2, a1, b1, c1, d1, a2, b2, c2, d2)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 50
   Unobserved stochastic nodes: 4
   Total graph size: 74

Initializing model

# 结果总结
summary(results)


Iterations = 5001:10000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

           Mean      SD Naive SE Time-series SE
meandiff 1.5962 0.37021 0.005236       0.005236
mu[1]    4.7853 0.21593 0.003054       0.003054
mu[2]    3.1891 0.30277 0.004282       0.004282
prob[1]  0.9998 0.01414 0.000200       0.000200
prob[2]  0.2656 0.44170 0.006247       0.006378
sdratio  0.8958 0.18977 0.002684       0.002847
tau[1]   0.7597 0.19559 0.002766       0.002803
tau[2]   0.5954 0.19581 0.002769       0.002915

2. Quantiles for each variable:

           2.5%    25%    50%    75% 97.5%
meandiff 0.8753 1.3456 1.5991 1.8463 2.314
mu[1]    4.3691 4.6403 4.7845 4.9260 5.209
mu[2]    2.5895 2.9936 3.1882 3.3913 3.790
prob[1]  1.0000 1.0000 1.0000 1.0000 1.000
prob[2]  0.0000 0.0000 0.0000 1.0000 1.000
sdratio  0.5697 0.7634 0.8789 1.0099 1.320
tau[1]   0.4174 0.6182 0.7424 0.8831 1.181
tau[2]   0.2737 0.4570 0.5724 0.7135 1.040

EXERCISE 5.30.

model_string <- "
model {
  for(i in 1:n[1]) {
    low[i] ~ dlnorm(mu[1], tau[1])
  }
  for(i in 1:n[2]) {
    normal[i] ~ dlnorm(mu[2], tau[2])
  }

  # 使用不同的先验分布
  mu[1] ~ dnorm(0, 0.00001) # 非信息性先验
  mu[2] ~ dnorm(0, 0.00001) # 非信息性先验
  
  tau[1] ~ dgamma(0.001, 0.001) # 非信息性先验
  tau[2] ~ dgamma(0.001, 0.001) # 非信息性先验
  
  med[1] <- exp(mu[1]) 
  med[2] <- exp(mu[2]) 
  rmed <- med[2] / med[1] 

  test[1] <- step(med[2] - med[1]) 
  test[2] <- step(Nf - Lf) 
  
  Lf ~ dlnorm(mu[1], tau[1]) 
  Nf ~ dlnorm(mu[2], tau[2]) 
  
  dmu <- mu[2] - mu[1] 
  rtau <- tau[2] / tau[1] 
}
"


data_list <- list(
  n = c(19, 15),
  low = c(91, 46, 95, 60, 33, 410, 105, 43, 189, 1097, 54, 178, 114, 137, 233, 101, 25, 70, 357),
  normal = c(370, 267, 99, 157, 75, 1281, 48, 298, 268, 62, 804, 430, 171, 694, 404),
  Lf = 50,
  Nf = 50
)


run_jags_model <- function(data_list, n.iter = 10000) {
  jags_model <- jags.model(textConnection(model_string), data = data_list, n.chains = 1, n.adapt = 5000)
  update(jags_model, n.iter = 5000)
  results <- coda.samples(jags_model, variable.names = c("mu", "tau", "med", "rmed", "test"), n.iter = n.iter)
  return(results)
}

# 运行模型
results <- run_jags_model(data_list)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 36
   Unobserved stochastic nodes: 4
   Total graph size: 54

Initializing model

# 结果总结
summary(results)


Iterations = 10001:20000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 10000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

            Mean      SD Naive SE Time-series SE
med[1]  108.7291 24.4553 0.244553       0.313422
med[2]  227.6707 65.1836 0.651836       0.833170
mu[1]     4.6644  0.2210 0.002210       0.002824
mu[2]     5.3893  0.2773 0.002773       0.003522
rmed      2.2013  0.8303 0.008303       0.010917
tau[1]    1.1430  0.3643 0.003643       0.003977
tau[2]    0.9707  0.3569 0.003569       0.003875
test[1]   0.9774  0.1486 0.001486       0.001849
test[2]   1.0000  0.0000 0.000000       0.000000

2. Quantiles for each variable:

            2.5%      25%      50%     75%   97.5%
med[1]   68.3571  91.9613 106.1596 122.213 165.776
med[2]  125.8918 183.1128 218.7977 262.320 381.281
mu[1]     4.2247   4.5214   4.6649   4.806   5.111
mu[2]     4.8354   5.2101   5.3881   5.570   5.944
rmed      1.0180   1.6316   2.0421   2.607   4.159
tau[1]    0.5271   0.8793   1.1068   1.377   1.938
tau[2]    0.4017   0.7141   0.9276   1.182   1.786
test[1]   1.0000   1.0000   1.0000   1.000   1.000
test[2]   1.0000   1.0000   1.0000   1.000   1.000

# 绘制结果（可选）
plot(results)

# 选择不同的先验分布
sensitivity_analysis <- function(mu1_prior, mu2_prior, tau1_prior, tau2_prior) {
  model_string <- sprintf("
  model {
    for(i in 1:n[1]) {
      low[i] ~ dlnorm(mu[1], tau[1])
    }
    for(i in 1:n[2]) {
      normal[i] ~ dlnorm(mu[2], tau[2])
    }

    mu[1] ~ %s
    mu[2] ~ %s

    tau[1] ~ %s
    tau[2] ~ %s
    
    med[1] <- exp(mu[1]) 
    med[2] <- exp(mu[2]) 
    rmed <- med[2] / med[1] 

    test[1] <- step(med[2] - med[1]) 
    test[2] <- step(Nf - Lf) 

    Lf ~ dlnorm(mu[1], tau[1]) 
    Nf ~ dlnorm(mu[2], tau[2]) 

    dmu <- mu[2] - mu[1] 
    rtau <- tau[2] / tau[1] 
  }
  ", mu1_prior, mu2_prior, tau1_prior, tau2_prior)

  jags_model <- jags.model(textConnection(model_string), data = data_list, n.chains = 1, n.adapt = 5000)
  update(jags_model, n.iter = 5000)
  results <- coda.samples(jags_model, variable.names = c("mu", "tau", "med", "rmed", "test"), n.iter = 10000)
  return(results)
}

# 运行敏感性分析示例
sensitivity_results <- sensitivity_analysis("dunif(-10000, 10000)", "dunif(-10000, 10000)", "dunif(0, 10000)", "dunif(0, 10000)")

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 36
   Unobserved stochastic nodes: 4
   Total graph size: 54

Initializing model

# 总结敏感性分析结果
 summary(sensitivity_results)


Iterations = 10001:20000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 10000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

            Mean      SD Naive SE Time-series SE
med[1]  108.3282 22.8235 0.228235       0.312057
med[2]  227.5207 59.1990 0.591990       0.804814
mu[1]     4.6638  0.2062 0.002062       0.002773
mu[2]     5.3952  0.2526 0.002526       0.003369
rmed      2.1934  0.7438 0.007438       0.010109
tau[1]    1.2723  0.3989 0.003989       0.005428
tau[2]    1.0999  0.3768 0.003768       0.005287
test[1]   0.9864  0.1158 0.001158       0.001464
test[2]   1.0000  0.0000 0.000000       0.000000

2. Quantiles for each variable:

            2.5%      25%     50%     75%   97.5%
med[1]   70.8756  92.3519 105.711 121.408 158.865
med[2]  133.3044 187.1056 219.895 259.969 362.284
mu[1]     4.2609   4.5256   4.661   4.799   5.068
mu[2]     4.8926   5.2317   5.393   5.561   5.892
rmed      1.0867   1.6707   2.075   2.592   3.963
tau[1]    0.6227   0.9863   1.231   1.506   2.160
tau[2]    0.4944   0.8320   1.053   1.326   1.964
test[1]   1.0000   1.0000   1.000   1.000   1.000
test[2]   1.0000   1.0000   1.000   1.000   1.000

plot(sensitivity_results)